
%%%%%%%%%%%%%%%%%%%%  Technical parameters  %%%%%%%%%%%%%%%%%%%%%%%%%
uptd=0.2; %Change to 0.1 if policy function fluctuate
outfreq=20; %Display frequency
iter_tol=600; % Iteration tolerance
tol=1e-5; %Numerical tolerance
%Number of States
NB=25;%21;%18; %number of nodes.
NB2=111;
NBG=25;%9;%7; %Number of nodes of government debt
bmin=.1; %-0.4; 
%bmin2=-.2; %-0.4; 
bmin2=-1.05; %-0.4;
bmax2=1.85;%1.75;
BGmin=0.0;
NT      = 6;                % Number of grid points for tradables
NR      = 1;                % Number of grid points for nontradables
%NKAPPA  =6;
NKAPPA  =5;
NPV  =5;
NNU  =1;
NSS=NR*NT*NKAPPA*NNU*NPV;

%CALIBRATED PARAMETERS
minNu=.0348;
maxNu=.1648;
kost0=(maxNu+minNu)/(2);
kost1=(maxNu-kost0)/.0475;
% CALIBRATED PARAMETERS
meank=0.476;
beta=0.954;
sigmak=0.0195;
v=0.0075;
p=0.91; % CORRELATION OF THE CHOICES UNDER REPAYMENT
theta=1/3; % PROBABILITY OF REENTRY

%CALIBRATED OUTSIDE
step_BG=0.05;
ita=(1/0.83)-1;
yn=1;
yN=ones(1,NSS);
std_kshocks=3;
std_PVshocks=2;
std_yshocks=3;
rbar=0.024895; %%1999-2012 German 1 Year bonds
spread_target=0.007251429; %Average 99-2012


euler_constant=0.5772156649; %EULER CONSTANT
mean_logi_shocks=-v*euler_constant;
sigma=2;
R    =  (rbar+1);
duration=6; %6years 
delta=(1+rbar+spread_target)/duration-rbar-spread_target; %Macayformula
%GENERATE INCOME SHOCKS TO TRADABLES

rho=.7534525;
sigmay=.010416;
[Z,Zprob] = tauchen(NT,0,rho,sigmay,std_yshocks); 
% SHOCKS TO PRIVATE DEFAULT

rho_spreads=.9135935  ;
sigma_spreads=.3242227  ;
pv_cte=-3.655621;
mean_spreads=exp(pv_cte);
[PVDefs,PVDefsprob] = tauchen(NPV,pv_cte,rho_spreads,sigma_spreads,std_PVshocks); 
%FINANCIAL SHOCK Kappa Shock to the Credit Constraint
rhok=rho;
[kappa_values,Prob_kappa]= tauchen(NKAPPA,meank,rhok,sigmak,std_kshocks);
Prob1=kron(PVDefsprob,Zprob);   % check this is the right
Prob=kron(Prob_kappa,Prob1);
yT_values      = exp(Z);
PV_values= exp(PVDefs);
PV_values0=kron(PV_values,ones(NT,NT))';
PV_values1=repmat(PV_values0,NKAPPA);
yT0=repmat(yT_values,NPV*NKAPPA)';
KAPPAS0=kron(kappa_values,ones(NT*NPV,NT*NPV))';
yT=yT0(1,:); %Exogenous Shocks Vector of Income
KAPPAS=KAPPAS0(1,:); %Exogenous shock Vector of Kappas
PVD=PV_values1(1,:); %Exogenous shock Vector of Default Cost
zT0=repmat(Z,NPV*NKAPPA)';
zT=zT0(1,:); 
NUS=max(zeros(1,NSS),kost0*ones(1,NSS)+ kost1*zT); %Exogenous shock Vector of Default Cost
clear zT0
clear yT0 KAPPAS0 PV_values0 PV_values1 Prob1 Z Zprob zT Prob_kappa PVDefsprob PVDefs mean_spreads
q=(ones(1,NSS)-PVD)*Prob'./R;
bmax_noG=(min(KAPPAS(:))+1)*min(yT(:))/(1-min(PVD(:)));
step_b=(1.3-.1)/(NB-1);
%% Calculate Omega
 [num,tex] = xlsread('targets.xls','Sheet1');
targets=[mean(num(:,5)), std(num(:,9)), mean(num(:,2)), mean(num(:,3)),std(num(:,6)), std(num(:,7)),mean(num(:,4)),std(num(:,8))];
ShareofT=.3938026;
Private_debt_target=targets(3);
Public_debt_target=targets(4);
Public_debt_factor=delta/(1-(1-delta)/R) ;
RatioOmega=mean((yT.*(1-ShareofT)).*(ShareofT.*(yT-Private_debt_target*rbar-delta*rbar*Public_debt_target/Public_debt_factor).^(1+ita)).^(-1));
omega=1/(1+RatioOmega);

%% GRID GENERATE
BGgrid=zeros(NBG,1);
BGgrid(1)=BGmin;
zeroBG_ind=1;

for i=2:NBG
    BGgrid(i)=BGgrid(i-1)+step_BG; 
end

bgrid=zeros(NBG,NB);

for x=1:NBG
    bgrid(x,1)=bmin;
end
for x=1:NBG
    for i=2:NB
        bgrid(x,i)=bgrid(x,i-1)+step_b; 
    end
end
bmax=max(max(bgrid));
BGmax=max(BGgrid);
shocksgrid=zeros(NSS,1);

for i=1:NSS
    shocksgrid(i)=i;
end
bmax_grid=max(max(bgrid));
bmin_grid=min(min(bgrid));

bgrid_large=zeros(NB2,1);

    bgrid_large(1)=bmin2;
bgrid_large=linspace(bmin2,bmax2,NB2);

for i=1:NBG
    for j=1:NB
        for t=1:NSS
          for inext=1:NBG  
            
            KAPPAS_big(i,j,t,inext)=KAPPAS(t);
        
             NUS_big(i,j,t,inext)=NUS(t);
             bp(i,j,t,inext)=bgrid(i,j);
             b_big(i,j,t,inext)=bgrid(i,j);
             BG_big(i,j,t,inext)=BGgrid(i);
             shocks_big(i,j,t,inext)=t;
             c_big(i,j,t,inext)=yT(t)+bgrid(i,j)*q(t)-(1-PVD(t))*bgrid(i,j)+(1/R)*(BGgrid(inext)-(1-delta)*BGgrid(i))-delta*BGgrid(i);
                BG_bignext(i,j,t,inext)=BGgrid(inext);
                yt_big(i,j,t,inext)=yT(t);
                   yn_big(i,j,t,inext)=yN(t);
                q_big(i,j,t,inext)=q(t);
                PVD_big(i,j,t,inext)=PVD(t);
         
          end
        end
    end
end


%% Matrices for The Sovereign Default Equilibruim 
%Policy function of Public Debt

V=zeros(NBG,NB,NSS);
Q=(1/R)*ones(NBG,NB,NSS);
V_upd=zeros(NBG,NB,NSS); 
PreQ=(1/R)*ones(NBG,NB,NSS);


qtemp_big=(1/R)*ones(NBG,NB,NSS,NBG);
evc_big=zeros(NBG,NB,NSS,NBG);
evdef_big=zeros(NBG,NB,NSS);
bgrid1d=squeeze(bgrid(zeroBG_ind,:));
bgrid2D=zeros(NB,NSS);
shocs2D=zeros(NB,NSS);
yT2D=zeros(NB,NSS);
nus2D=zeros(NB,NSS);
yN2D=zeros(NB,NSS);
q2D=zeros(NB,NSS);
PVD2D=zeros(NB,NSS);
for i=1:NB
      for j=1:NSS
      shocs2D(i,j)=j;
      yT2D(i,j)=yT(j);
      nus2D(i,j)=NUS(j);
      q2D(i,j)=q(j);
      PVD2D(i,j)=PVD(j);
      yN2D(i,j)=yN(j);
      bgrid2D(i,j)=bgrid1d(i);
         KAPPAS2D(i,j)=KAPPAS(j);
      end
end

    for i=1:NB
             PVDm1_2d(i,:)=1.-squeeze(PVD);
          end